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Communicated  by  E.  Y.  Rodin 


Abstract — Based  on  the  Newton-rRaphson  method,  this  paper  presents  recursive  algorithms  that  are 
rapidly  convergent  and  more  stable  for  modeling  the  equivalent  continuous-time  (discrete-time)  model 
from  the  available  discrete-time  (continuous-time)  model  for  a  fixed  sampling  period.  The  newly  developed 
recursive  algorithms  relax  the  constraints  imposed  upon  the  existing  model  conversion  algorithms,  and. 
thus,  enhance  the  applications  of  microprocessors  and  associated  microelectronics  to  digital  control 
systems.  A  practical  example  is  presented  to  demonstrate  the  effectiveness  of  the  proposed  procedures. 


5. 


1.  INTRODUCTION 


A  practical  dynamic  system  is  often  described  by  a  sampled-data  system  which  is  a  collection  of 
continuous-time  and  discrete-time  models  with  their  own  characteristics.  For  simplicity  in  analysis 
and  design  the  continuous-time  (discrete-time)  subsystems  are  converted  into  equivalent  discrete¬ 
time  (continuous-time)  subsystems  so  that  the  original  sampled-data  system  can  be  completely 
analyzed  and  designed  in  the  discrete-time  (continuous-time)  domain  via  the  well-developed 
techniques.  Hence,  the  conversion  of  the  continuous-time  (discrete-time)  model  to  the  equivalent 
discrete-time  (continuous-time)  model  is  often  required.  Advances  in  industrial  electronics  and 
computer  technologies  have  made  a  dramatic  extension  of  the  possibilities  of  control  algorithms 
to  be  implemented  via  these  high  performances,  low  cost  microprocessors  and  associated 
microelectronics;  however,  microprocessors  are  slow  digital  machines,  and  are  usually  equipped 
with  a  small  wordlength.  For  effective  use  of  microprocessors  in  digital  control  systems,  the 
selection  of  a  suitable  sampling  time  and  the  development  of  recursive  algorithms,  which  possess 
fast  convergence  rate  and  numerical  stability  to  offset  the  slow  computational  speed  and  round-off 
errors,  become  an  important  part  of  any  digital  control  system.  The  choice  of  the  sampling  time 
for  the  identification  and  model  conversion  of  control  systems  has  been  discussed  in  Refs  [1,2], 
For  example,  a  rule  of  thumb  for  identifying  a  continuous-time  system  from  samples  of 
input/output  data  has  been  suggested  in  Ref.  [3],  which  states  that  the  sampling  interval  T  should 
be  chosen  in  such  a  way  that  |  A„r|  <  0.5,  where  X„  is  the  magnitude  of  the  largest  eigenvalue  of 
the  continuous-time  model.  Without  knowing  the  largest  eigenvalue  of  the  continuous-time  system 
to  be  identified,  a  simple  procedure  for  obtaining  a  good  choice  of  T  has  been  proposed  in  Ref. 
[4].  For  critical  dynamic  systems,  such  as  the  power  systems  [5]  and  missile  systems  (6,  7]  which 
contain  large  and  small  eigenvalues,  a  very  small  sampling  period  T  is  usually  required  for  the 
identification  and  control  of  the  systems.  However,  for  the  use  of  the  identified  model  in  digital 
control  systems,  the  sampling  period  T  cannot  be  chosen  as  small  as  we  would  like  due  to  the 
consideration  of  the  execution  time  of  a  particular  microprocessor  and  the  round-off  errors  caused 
by  a  small  wordlength  of  the  microprocessors.  In  other  words,  the  choice  of  the  sampling  period 
T  depends  upon  not  only  the  dynamic  property  of  each  subsystem  but  also  the  constraints  of  the 
microprocessors  used  for  digital  control  systems.  In  all,  the  development  of  a  rapidly  convergent 
and  numerically  stable  recursive  algorithm  with  a  fixed  sampling  period  for  the  model  conversions 
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is  necessary  in  order  to  enhance  the  applications  of  microprocessors  and  associated  micro¬ 
electronics  to  digital  control  systems. 

Based  on  the  Newton-Raphson  method,  this  paper  develops  rapidly  convergent  and  more  stable 
recursive  algorithms  with  a  fixed  sampling  period  for  continuous-time  and  discrete-time  model 
conversions.  The  proposed  algorithms  enable  us  to  reduce  the  constraints  imposed  upon  the 
existing  model  conversion  techniques  and.  thus,  enhance  the  applications  of  microprocessors  and 
associated  microelectronics  to  digital  control  systems.  In  Section  2.  we  review  some  non-recursive 
algorithms  for  finding  the  equivalent  continuous-time  (discrete-time)  model  from  an  available 
discrete-time  (continuous-time)  system.  In  Section  3,  we  develop  a  rapidly  convergent  and 
numerically  stable  recursive  algorithm  with  a  fixed  sampling  period  T  for  modelling  continuous¬ 
time  model  from  a  discrete-time  model.  In  Section  4.  we  also  develop  a  rapidly  convergent  and 
numerically  stable  recursive  algorithm  with  a  fixed  sampling  period  for  finding  the  discrete-time 
model  from  a  continuous-time  system.  A  practical  example  is  given  in  Section  5  and  the  results 
are  summarized  in  Section  6. 

2.  REVIEW  OF  SOME  NON-RECURSIVE  ALGORITHMS  FOR 
MODEL  CONVERSIONS 

There  are  many  model  conversion  techniques  [1,  2,  8-12]  available  in  the  literature.  The  highly 
recommended  techniques  are  summarized  as  follows. 

{a)  The  scaling,  squaring,  continued-fraction  method  for  finding  a  discrete -time  model  from  a 
continuous -time  model  [8,  10] 

Let  the  continuous-time  system  be 

X{t)  =  AXit)-i-Bu(t);  ;)r(0),  (1) 

where  Xf)  and  uf)  are  the  n  x  1  state  vector  and  the  m  x  1  input  vector,  respectively,  and  the 
A  and  B  are  constant  matrices  of  appropriate  dimensions.  If  we  approximate  «(t)  as  a  piecewise 
input  function 

u(t)  =  u{kT) -- u(k)  for  kT^t<{k-h\)T  (2) 

where  the  T  is  the  sampling  period  with  T  <  njto,  where  the  to,  is  the  highest  frequency  component 
in  equation  (1),  then  we  can  write  the  equivalent  discrete-time  model  as 

2r(A:  -F  1)  =  GX{k)  -f  Hu{ky,  2r(0),  (3a) 

where 

G  =  exp(/4r)  (3b) 

and 

H  =  IG  -f]A-'B.  (3c) 

The  matrix  f„  denotes  an  n  x  n  identity  matrix  and  the  vector  X{k)  =  X(kr)  =  X(t)  t  —  kT.  The 
constraint  of  the  sampling  period  r(<7t/aj,)  in  equation  (2)  implies  that  the  obtained  discrete-time 
system  matrix  G  is  a  nonsingular  matrix  without  any  negative  real  eigenvalues.  The  matrices  G  and 
H  can  be  determined  exactly  from  the  matrices  A  and  B  using  the  eigenvalues  and  eigenvectors 
approaches  [8].  However,  for  simplicity,  in  computations  without  explicitly  involving  the  eigen¬ 
values  and  eigenvectors  of  the  matrix  A  which  may  be  defective,  the  approximant  of  the 
matrix-valued  function  in  equation  in  (3b),  exp(>4r),  is  obtained  for  modelling  the  matrix  G.  The 
most  commonly  used  method  for  determining  the  approximant  is  the  direct  truncation  of  the 
following  infinite  series: 

G=tx^(AT)  =  f  +  AT  +  \-MTf^y  =  Y^\-{ATy.  (4) 

The  truncation  of  the  above  series  results  in  a  good  approximant  if  r«  1  and  A  is  not  a  stiff  matrix. 
A  stiff  matrix  is  defined  as  a  matrix  which  contains  both  large  and  small  eigenvalues.  The  analysis 
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of  the  truncation  error  for  the  series  approximation  in  equation  (4)  has  been  discussed  in  Ref.  [13], 
while  the  round-off  errors  have  been  studied  in  Ref.  [14].  For  matching  the  sampling  period  of  the 
existing  discrete-time  subsystems  in  the  original  sampled-data  system  so  that  the  modified 
sampled-data  system  has  a  uniformly  sampled  rate,  the  sampling  period  T  in  equation  (4)  can  not 
be  arbitrarily  chosen.  For  a  fixed  sampling  period,  a  popular  method  for  determining  the  matrix 
G  from  exp(^r)  is  the  scaling,  squaring  continued-fraction  method  [8]  as  follows: 

G  =  (expiATIq))" 

^[f  +  ATIqY  ^G,,, 

{[f-{(AT/q)]-‘[f  +  -2iATIq)]y^  Gy, 


(5a) 

(5b) 

(5c) 


In  general,  the  approximants  of  exp(^D  can  be  written  as 


G,,,=  ^ 


fPV^'  (/-l-y)!(//2-l)!  / 

ATV'- 

lLv  =  o(i-l)!7!(<72-l-7)!\ 

9  j. 

for  i  =  2, 4,  6, . . . 

(n.-^v2  (,  _i  _;)!((/ _n/2)!  / 

'  ATy 

iL/.o  (i-l)!7!((i-  l)/2-;)!\ 

.  9  j 

for  I  =  3,  5,  7, . . . 

(/- l-y)!(//2)!  /ATV 


L?o(<-i)!yU»72-i-;)!V  9 


(5d) 


(/_i  _;)!((,•  _i)/2)!  /AT\ 


Ijto  (/-l)!y!((» -l)/2-y)!V  9  J  1 


(5e) 


The  first  subscript  i  of  G,.,  in  equation  (5)  denotes  the  number  of  the  quotients  in  the  continued 
fraction  expansion  of  exp(^r/<?)  to  be  taken  and  the  second  subscript  ^  of  G,.,  is  the  scaling  and 
squaring  factor  which  is  often  chosen  as  a  power  of  two.  For  the  analysis  of  the  truncation  error 
and  the  determination  of  the  q,  we  define  a  relative  error  matrix  and  compute  its  approximant  as 


where 


Eo^[G,,-G]G-' 

-  9  —  ) 

\<}  Jl 

(6a) 

r(_iy/2,!(/_  1)! 

for  i  =  2,4,6, ... , 

(6b) 

(//2)!(//2  -  1)!  ’ 

,  for  i  =  3,5,1, ... . 

(6c) 

L((/-l)/2)!((/-l)/2)!^ 

With  the  prescribed  values  such  as  the  sampling  period  T,  the  value  i  in  equations  (5)  and  (6)  and 
the  relative  error  tolerance  ll^dl  in  equation  (6a),  we  can  determine  the  scaling  and  squaring  factor 
q  in  equations  (5)  from  (6)  as 

^  =(ll(/iryii/(|/i:|i|£cll))''‘-”' 


Note  that  a  small  l|£cll  results  in  a  large  q.  When  the  system  matrix  A  in  equations  (5)  is  a  stiffs 
matrix  and  the  power  j  is  a  large  value,  the  computation  the  high  order  approximant  in^ 
equations  (5)  often  produces  round-off  errors  due  to  the  cc  .’Du'ation  of  [ATjqy.  To  avoid  the 
above  numerical  difficulty,  we  shall  develop  a  stable  and  fast  ’rsive  algorithm  for  computing 
the  matrix  G  from  the  matrix  A  in  Section  4. 


□ 

□ 


(b)  The  square-root,  scaling,  continued-fraction  method  for  finding  a  continuous-time  model  from  a  _ 
discrete -time  model  [12] 

The  equivalent  continuous-time  model  in  equation  (1)  can  be  obtained  from  the  discrete-time  ^ 
model  in  equations  (3)  as  follows: 


A  =-ln(G) 


7  CodfiS 


(8a)  nd/or 
dl 


u 


f('\  ^ 


o 
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and 


B  =  A[G (8b) 

The  equivalent  continuous-time  system  matrix  A  in  equation  (8a)  can  be  obtained  by  truncating 
the  following  infinite  series  as 


—  y  In  (G)  —  j  -I-  5  -1- 1  /?’-(- 


=  -  Y  ~ 


-H  1 


R- 


(9) 


where  R  =  [G  —  I„\{G  4-/,]''  with  Re[<T(G)]>0,  where  a{G)  is  the  spectrum  of  the  matrix  G. 

When  Re[<T(G)]  ->  -1-  1  and  Im[<T(G)]  -►0.  then  |(t(/?)|«1.  Hence,  the  truncation  of  the  above 
series  will  give  a  good  approximant.  Since  not  all  eigenvalues  of  the  matrix  G  always  lie  in  the  right 
half  complex  plane  (i.e.  Re[<T(G)]  >  0)  and  close  to  a  positive  unity;  therefore,  we  need  to  modify 
the  matrix  G  before  we  carry  out  the  truncation  of  the  infinite  series  in  equation  (9).  Taking  the 
(^th  root  of  the  matrix  G  in  equation  (5a)  results  in 

G  =  [exp(.4  r)]'"'  =  exp(/l  Tjq).  (10a) 


Thus, 


A 


ln(G'''') 


with  ^  as  a  power  of  two. 


(10b) 


If  the  matrix  G''"  is  the  principal  ^th  root  [15,  16]  of  the  non-singular  matrix  G,  where  the  principal 
^th  root  of  a  matrix  G  has  the  properties  that  (G'  ’)’  =  G  and  arg((T(G ''"))£(  —  -i-t:/^)  for  <7  >2, 

then  the  matrix  A  in  equation  (10b)  can  be  approximated  by  truncating  the  following  infinite  series 
as  follows: 


A=]j.\n{G)  =  ^\n(G'’%  (11a) 

where 

ln(G''’)  =  2[R-fiR'-(-iR'  +  ---]  =  2  Y 

j.o 2;  +  1 

=:2R=2[R-f]R’]s:  •  •• 
and 

For  a  large  q,  Re[<T(G''’)]  -►  -I-  1  and  Im[ff(G '''')]  ->  0.  Since  q  is  chosen  as  the  power  of  two,  the 
principal  ^th  root  of  the  matrix  G  can  be  found  by  successively  taking  its  principal  square  root 
of  the  matrix  G.  A  numerically  stable  recursive  algorithm  [17,  18]  with  a  quadratic  convergence  rate 
can  be  applied  to  determine  G'"  and  is  summarized  as  follows: 


F(/c-f-l)  =  nF(k)-he(A:)-'], 

P{0)  =  G 

(12a) 

C(*  +  l)  =  HF(k)-'  +  C(A:)], 

2(0)  =  /„ 

(12b) 

lim  P{k)  =  G''^ 

*  -*  X 

(I2c) 

lim  Q{k)  =  (G''Y'- 

(12d) 

k  —  00 


The  algorithm  in  equations  (12)  fails  [16,  17,  18]  when  the  matrix  G  contains  any  negative  real 
eigenvalues.  In  this  case  the  matrix  G  can  be  rotated  by  a  small  positive  real  angle  Ad  to  give 
G  =  G  exp(-yA0)  so  that  G’'^  =  G''^exp(yA0/2). 

From  equation  (10a)  we  observe  that  when  the  continuous-time  system  matrix  AT  is  normalized 
by  a  scaling  factor  q,  the  equivalent  discrete-time  matrix  becomes  G Therefore,  the  scaling  factor 
q  for  the  continuous-time  system  matrix  becomes  the  square -root  factor  for  the  discrete-time  system 


(Hb) 

(11c) 
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matrix.  Also,  from  equation  (10b)  it  can  be  seen  that  the  scaling  factor  q  for  the  discrete-time  system 
matrix  becomes  the  squaring  factor  for  the  continuous-time  system  matrix  in  equation  (5a).  Once 
we  have  the  matrix  G' ",  we  can  determine  the  approximant  of  ln(G' ")  via  the  continued-fraction 
method  [12].  Some  approximants  of  InfC'")  obtained  by  the  square-root,  scaling,  continued- 
fraction  method  can  be  written  as  follows; 

>1  =yln(G)  =  |;ln(G'’),  (13a) 

where 

(13b) 

^j2R[In-\Rr'  (13c) 

a  a,,  (13d) 

where  the  matrix  R  is  defined  in  expression  (lid). 

The  first  subscript  /  of  denotes  that  i  quotients  of  the  continued-fraction  expansion  of  ln(G ' ") 
have  been  used  for  finding  the  approximant  of  InfC'"),  whereas  the  second  subscript  q  of  A,  ^  is 
the  square-root  and  scaling  factor.  For  the  analysis  of  the  truncation  error  and  the  determination 
of  the  q,  we  define  a  relative  error  matrix  and  compute  its  approximant  as 

Ra  -  [^A,„  -A]A-'c=,  -£,^2'  ~  L[G  -f  G-'  -  lIJKlqf,  (1 4a) 

where 

rni-\)i2 

n  [(/  +  2)\i  +4y...  (2i  -  1)^(2/  +  1)]-',  for  /•  =  1,  3,  5 . 

L  /-o 

L^<  (14b) 

I  fl  +  !)=(/  +  3)^  . .  (2i  -  1)^(2/  +  1)]- for  i  -  2,  4,  6 . 

'-Lz-O  - 

Having  prescribed  the  values  for  the  sampling  period  T,  the  value  of  i  in  equation  (13)  and  the 
relative  error  tolerance  ||£^||  in  equations  (14),  the  square-root  and  scaling  factor  q  can  be 
determined  as  follows: 

q  +G-'  -2I„y\\l\\EA\)''^-.  (15) 

A  small  llE^II  results  in  a  large  q.  When  the  matrix  G  is  a  stiff  matrix,  the  computation  of  G'" 
with  a  high  power  q  in  equations  (12)  may  induce  round-off  errors.  To  overcome  the  above 
difficulties,  the  development  of  a  fast  and  stable  recursive  algorithm  for  finding  the  matrix  A  from 
the  matrix  G  becomes  necessary  and  is  shown  as  follows. 

3.  FAST  AND  STABLE  RECURSIVE  ALGORITHM  FOR  FINDING  A  FROM  C 

The  recursive  algorithm  developed  in  Refs  [1, 2]  for  finding  the  matrix  A  from  the  matrix  G  can 
be  reviewed  as  follows. 

The  desired  matrix  A  can  be  solved  from  the  following  matrix  equation; 

F,(/I)  =  exp(-/ir)-G-' =0„  (16a) 

where  the  matrix  0„  is  an  n  x  n  null  matrix.  Following  the  Newton-Raphson  method  yields 

Aik +  \)  =  A{k)-  [F'iA ik ))] -'FiAik)).  (1 6b) 

Substituting  the  Jacobian  matrix,  F’iA{ky)  =  —T  cxpi  —  Aik)T)  which  is  a  non-linear  matrix- 
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valued  function,  in  equation  (1 6b)  gives 

/l(fc  +  l)  =  .4(/c)  +  -^[/,-exp(/t(A:)r)G-'],  for  A:  =  0.  1,  2, . .  . ,  (16c) 

^(0)  =  ^2(G-/J(G +  /,)-’  (16d) 

lim  /l(/c)  =  A.  (16e) 

i  —  X 

At  each  iteration,  the  approximant  of  exp{A(k)T)  via  the  direct  truncation  method  in  equation 
(4)  is  determined  and  substituted  into  equation  (16c)  for  finding  the  up-dated  y4(A:  -t-  1).  In  Refs 
[1,2,  11],  it  was  stated  that  the  desired  matrix  A  can  be  determined  when  |<t(/1  (/c))ri  ^0.5  and 
the  algorithm  in  equations  (16)  may  become  unstable  for  a  given  matrix  G  if  the  sampling  period 
T  is  too  large.  The  inverse  of  the  matrix  G  in  equation  (16c)  can  be  avoided  if  we  use  an  alternative 


matrix  equation  given  as 

^2(^)  =  — exp(/ir) -f  G  =  0„.  (17a) 

Thus,  the  corresponding  recursive  algorithm  becomes 

^(A: -I- l)  =  ^(A:)-i[/„-exp(-,4(A:)r)G],  for  A:  =0,1,2 .  (17b) 

^(0)  =  |(G-/,)(G +/„)-'  (17c) 

lim  /4(A:)  =  /<.  (17d) 

k-  X, 


The  accuracy  of  the  Newton-Raphson  algorithms  shown  in  equations  (16)  and  (17)  is  heavily 
dependent  on  the  approximant  of  exp(/l(A:)r)  or  exp(  — /4(A:)r).  In  this  section,  we  shall  develop 
a  rapidly  convergent  and  numerically  stable  recursive  algorithm  for  the  model  conversion  of  the 
system  which  may  not  satisfy  the  conditions  Re[a(G)]>0  and/or  lCT(/l)ri  <  0.5  with  a  fixed 
sampling  period  T. 

It  is  well  known  that  the  convergence  of  the  Newton-Raphson  algorithm  depends  heavily  upon 
the  initial  guess  of  the  recursive  algorithm,  and  the  property  of  the  Jacobian  matrix  (the  slope)  of 
a  matrix-valued  function  (a  scalar-valued  function).  For  example,  we  consider  two  scalar-valued 
non-linear  functions  f^{a)  =  e\p(  —  aT)~g~'  and  f2(a)=—exp{aT)  +  g  where  T  and  g  are 
constants  and  a  is  a  variable.  The  functions  /i(a)  and  /2(a)  have  the  slopes  —T/g  and  —Tg, 
respectively,  at  the  solution  a  =  l/T\n{g).  In  this  paper,  we  propose  a  new  function 
h{a)  =  ln(exp(  — ar)g)  =  -aT  -t-  In  (g)  which  is  a  linear  function  with  a  constant  slope  —  T  for  any 
a  including  the  solution  a  =  l/rin(g).  To  compare  the  convergence  speeds  of  the 
Newton-Raphson  algorithms  developed  for  the  non-linear  functions  /,(«)  and  fiia)  with  the 
convergence  speed  for  the  linear  function  h{a),  we  modify  the  functions  /i(a)  and  /2(a)  as 

/i(a)  =  (exp(-ar)-g'')g  =exp(-ar)g  -  1, 
and 

fiia)  =  +  g)g~'  =  -exp(ar)g-'  +  1, 

respectively,  so  that  the  modified  non-linear  functions  /i(a)  and  /2(a)  have  the  same  slope,  —  T, 
as  that  of  h(a)  at  the  solution  a  =  l/Tlnlg).  The  Inunctions  /,(fl)  and  /2(a)  are  indeed  the 
approximants  of  the  newly  proposed  function  h(a)  as  follows: 

/i(a)  =  ln(exp(-a7’)g)  =  ln((exp(-ar)g  -  1)+  1] 

=  [exp(-ar)g  -  l]-nexp(-ar)g  -  1]^ 

+  Hexp(-flr)g  -  Ip - 2:exp(-ar)g  -  1  =/i(a). 


(18a) 
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Also,  in  the  same  fashion,  we  obtain 

/i (a )  =  ln[exp(- a T)g]  =  -ln[exp(a7)^'']  ~  -exp(ar)g-'  +  1  =A(a).  (18b) 

Another  approximant  of  the  function  h{a),  defined  as  /i|(a),  can  be  obtained  as  follows; 

/i(a)  =  ln[exp(-ar)g]  =  2{[(exp(-ar)g  -  l)(exp(-ar)g  +  1)-'] 

+  jKexpf-aHg  -  l)(exp(-a7')g  +  !)]’  +  •  ■  •} 

=;  2(exp(-a7')g  -  l)(exp(-a7')g  +  1)-'  A  (I8c) 

In  addition,  an  alternative  approximant  of  the  function  h{a),  defined  as  /ii(a),  can  be  written  as 
h{a)  =  ^[expf-aDg] 

=  ln[(exp(-ar)g)''V(exp(-ar)g)“'  =:  (exp(-ar)g)''- 

-  (exp(-a7’)g)-''’ A (I8d) 

Note  that  all  functions  in  equations  (18)  have  the  same  solution  at  a  =  \/T\n(g)  and  their 
respective  slopes  for  any  a  are 

/;(a)=  -rexp(-ar)g, 

/2(a)  =  -Texp(aT)g'\ 

h'(a)=-T, 

h\(a)  =  -AT  txp{-aT)g{txp{-aT)g  +  1)"^ 

and 

A^a)  =  -  r[(exp(-ar)g)''^  +  (exp(-ar)-''^/2. 

The  curves  of  /i(a),  /2(a),  h{a),  h^{a)  and  /12(a)  are  shown  in  Fig.  1.  From  Fig.  1,  we  observe  that 
the  curve  h{a)  is  a  straight  line  and  the  curves  A, (a)  and  hi(a)  are  closer  to  the  curve  A  (a)  than 
those  of  /i(a)  and  /2(a).  Thus,  if  we  can  determine  the  exact  values  of  A  (a),  A,  (a)  and  /(a)  for 
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(■  =  1,  2  at  an  arbitrary  initial  guess  a  =  floi  then,  we  shall  obtain  the  solution  a  =  XjT  ln(g)  at  one 
step  for  the  linear  function  h{a),  a  few  steps  for  the  functions  /i,(a)  and  many  steps  for  the  functions 
/,(a)  with  i  =  1,2.  Thus,  the  convergence  speeds  of  the  Newton-Raphson  algorithms  developed 
for  the  functions  h,(a)  are  faster  than  those  of  the  functions  /.(a)  for  /  =  1, 2.  In  other  words,  the 
convergence  speed  of  the  Newton-Raphson  algorithm  depends  on  the  approximation  of  the 
logarithmic  function,  ln[exp(  — uTIg],  and  the  accuracy  of  the  recursive  algorithm  depends  on  the 
approximant  of  the  exponential  fui.ction,  exp(  — aT). 

Based  on  the  above  reasoning,  we  propose  a  new  matrix-valued  function  H{A)  to  improve  the 


recursive  algorithms  shown  in  equations  (16)  and  (17)  as  follows: 

/f(^)  =  ln[exp(//r)(;-']=  -ln[exp(-/4r)G]=  -ln[/J  =  0„.  (19a) 

The  corresponding  Newton-Raphson  algorithm  becomes 

>l(fc  +  1)  =  (19b) 

=  /^(^)-^ln(exp(/«(A:)r)(7-'],  for  0,1,2 _  (19c) 

=  A(k)  +  ]^\n[exp{-A(k)T)Gl  for  jt  =  0,  1,  2, . . . ,  (l9d) 

^(0)  =  0„  (19e) 

\{mA{k)  =  A.  (190 

k-co 


Note  that  the  first  iteration  yields  /4(1)  =  l/rin(G).  Because  the  Jacobian  matrix,  H'(A(k))  =  TI„, 
is  a  constant  diagonal  matrix,  the  convergence  speed  of  the  recursive  algorithm  in  equation  (19d) 
depends  on  the  approximant  of  ln[exp(  — /l(/:)r)G],  whereas  the  accuracy  of  the  recursive 
algorithm  depends  on  the  approximant  of  t%p(  —  A(k)T).  Since 

ln[exp(  — /( (A:)r)G]  ^  I„  —  exp(/l(fc)r)G"', 

the  recursive  algorithms  in  equations  (16)  and  (17)  are  a  special  case  of  the  algorithms  in  equations 
(I9d)  and  (1 9c),  respectively. 

The  logarithmic  matrix-valued  functions  in  equations  (19c)  and  (19d)  can  be  approximated  via 
either  the  direct  truncation  method  in  equation  (11)  or  the  square-root,  scaling  and  continued- 
fraction  method  in  equation  (13)  with  q'^\.  For  example,  when  q  =2  and  the  simplest 
approximant  in  equations  (11)  or  (13)  is  used,  we  obtain 

A(k +  \)  =  A{k)  +  ^[txp(-A(k)T)G  -  /,] [exp( - A(k)T)G  +  Q-' ,  (0)  =  0„  (20) 

where  A(k)^  A(k)/2  and  G  ^  C'^.  The  G  can  be  determined  via  the  stable  algorithm  in  equations 
(12).  The  direct  truncation  method  in  equation  (4)  or  the  scaling,  squaring  and  continued-fraction 
method  in  equation  (5)  can  be  applied  to  determine  the  approximant  of  txp(  —  A(k)T).  Note  that 
for  the  algorithm  in  equation  (20),  the  first  iteration  gives  /4(l)  =  2/7" [G  —  I„]{G  -t-  which  is 
the  initial  condition  used  in  the  algorithms  in  equations  (16)  and  (17).  It  can  be  shown  via  the 
practical  example  in  Section  5  that  both  algorithms  in  equations  (16)  and  (20)  are  numerically 
unstable.  To  overcome  the  instability,  a  new  recursive  algorithm  is  developed  for  finding  A  from 
G  as  follows: 

Define  M(k)  ^  cxp(  — A(k)T)G.  From  equation  (19d)  we  can  derive  the  equations  as  follows: 
exp(-i(fc  -I-  DDG  =exp(-^(k)7')G  exp(-ln(exp(-i(k)r)G)), 
i.e. 

M(k +  0  =  M{k)ixp(-\n(M(k))),  A/(0)  =  G''‘'  (21a) 

and 

i(k  +  l)  =  /i(A:)-|-lln(A/(A:)),  i(0)  =  0„  (21b) 
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\\mM{k)  =  I„  (21c) 

A  -•  X 

lim  qA{k)  =  A.  (21d) 

k  -*  x. 

The  following  approximants  can  be  used  to  approximate  \n(M(k))  in  equation  (21): 

\n{M{k))T-.M(k)~I„  (22a) 

~[M{k)-l„]-\[M^k)-I„\'  (22b) 

» 

or 

ln(A/(ilc))=r2/?(A:)  (22c) 

^  2R(k)[I„  -  -  \R{k)-]-'  (22d) 


where 

R(k)-^mfc)-/„]mk)-i- /„]-'■ 

If  the  approximant  in  expression  (22a)  is  used,  the  algorithm  in  equations  (21)  becomes 


M(ik  +  1)  =  A/(ifc)exp[/,-M(it)],  Af(0)  =  G''‘'  (23a) 

A(,k  +  \)  =  A{k)~j[I„-M{k)l  ^(0)  =  0„  (23b) 

\\mM(k)  =  J,  (23c) 

*  -•  00 

[imqAik)  =  A.  (23d) 

k  —  X 

If  the  approximant  in  expression  (22b)  is  used,  the  algorithm  in  equations  (21)  becomes 

yV/(/c  +  l)  =  M(*)exp[(/,-M(ik))  +  i(/„-A/(/c))^,  M(,0)  =  G'^^  (24a) 

A(k  +  l)  =  A(k)-^l(I„-M{k))  +  'i(I„-Mik)n  A(0)  =  0„  (24b) 

lim  M(k)  =  I„  (24c) 

fc  -•  X 

lim  ^/4(/c)  =  ^.  (24d) 

*-*« 

If  the  approximant  in  equation  (22c)  is  utilized,  the  algorithm  in  equations  (21)  becomes 

M(k  +  I)  =  A/(/fe)exp[-2/?(/t)],  M(0)  =  G'>‘>  (25a) 

A{k  +  \)  =  A{k)  +  ~R(k),  ^'(0)  =  0„  (25b) 

lim  M{k)  =  I„  (25c) 

k  ^  <x} 

lim  qA(k)  =  A.  (25d) 

k  —  aO 

Also,  if  the  approximant  in  expression  (22d)  is  employed,  the  algorithm  in  equations  (21)  gives 
A/(A:  +  l)  =  A/(*)exp{-2/?(*)[/„-i^/?(^)^[;„-j/?(*)^-'},  A/(0)  =  G''’  (26a) 

A(k  +  l)  =  A{k)  +  jR(k)[I„-fiR{kmh-\R{ky]-\  i(0)  =  0„  (26b) 

lim  M{k)  =  I„  (26c) 

k  -*  00 

lim  qA{k)  =  A.  (26d) 

k  ^  ao 
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The  matrix  exponential  functions  in  equations  (23a).  (24a).  (25a)  and  (26a)  can  be  computed  via 
the  techniques  shown  in  Section  2,  while  the  initial  conditions  C''  can  be  obtained  using  the 
algorithm  in  equations  (12).  The  algorithms  in  equations  (23).  (24),  (25)  and  (26)  are  numerically 
stable  and  the  stability  is  proved  in  Appendix  A.  Note  that  the  accuracy  of  the  recursive  algorithm 
in  equations  (21)  depends  upon  the  approximation  of  the  matrix  exponential  function,  while  the 
convergence  speed  of  the  recursive  algorithm  in  equations  (21)  depends  on  the  approximation  of 
the  matrix  logarithmic  function  and  the  factor  q. 

Substituting  the  required  system  matrix  A  and  the  available  matrix  G  to  equation  (8b)  yields  the 
input  matrix  B.  Thus,  we  obtain  the  desired  continuous-time  model  in  equation  (1). 


4.  FAST  AND  STABLE  RECURSIVE  ALGORITHM  FOR  FINDING  G  FROM  A 
From  the  newly  developed  function  in  equation  (19a)  for  finding  A  from  G.  we  can  construct 


the  following  matrix  equation  for  finding  G  from  A : 

E(G)  =  exp(]n[exp(-AT)G]}  -exp[ln(/J]  =  exp[ln(G)-.-<r]  -  /,  =  0,.  (27a) 

The  corresponding  Newton-Raphson  algorithm  becomes 

G(k  +  l)  =  G(/c)-[£-(G))-'£(G)lc.«*,  (27b) 

=  G(A:)exp[/ir-ln(G(/t))],  for  A:  =0.1,2 .  (27c) 

G(0)  =  /„  (27d) 


and 

Urn  Gik)  =  G.  (27e) 

k  X 

Note  that  the  first  iteration  gives  G(l)  =  exp(v*ir).  The  Jacobian  matrix, 

E'(G)  =  {exp[ln(G)  -  AT]]G''  exp(  -AT) 

is  a  constant  matrix.  Therefore,  the  co'^vergence  speed  depends  upon  the  approximant  of 
exp(  — /IT"  +  ln(G(/c))]  and  the  accuracy  of  the  recursive  algorithm  depends  on  the  approximant  of 
ln(G(Ar)).  Because  the  series  expansion  of  ln(G)  in  equation  (9)  for  Re[(7(G)]>0  converges  very 
fast,  good  approximants  of  ln(G(A:))  can  be  determined  even  if  the  matrix  G  is  a  stiff  matrix.  As 
a  result,  a(AT  —  ln(G(/c))]  0  and  good  approximants  of  exp[,4r  —  ln(G(A:))l  can  be  obtained  via 

the  methods  shown  in  equation  (4)  and  (5).  The  imaginary  parts  of  the  eigenvalues  oi  AT  may  be 
such  that  n  >  \\m(a(AT))  \  ^  njl.  As  a  result,  Re((T(G))  may  be  negative  and  the  infinite  series 
expansion  of  ln(G)  in  equation  (9)  will  not  converge.  To  overcome  this  difficulty,  we  normalize  the 
matrix  /I  by  a  scaling  factor  q  ^  Iso  that  nil  >  |Im((T(/l  T))!  >  0  and  Refer  (G))  >  0,  where  A  ^  A  Iq 
and  G  —  expf/^T).  Thus,  the  infinite  series  expansion  of  ln(C)  in  equation  (lib)  converges.  The 
Newton-Raphson  algorithm  in  equation  (27c)  becomes 

G(A- -F  1)  =  G(A)exp[/ir- ln(G(A:))],  for  A  =0,1,2 .  (28a) 

G(0)  =  /,  (28b) 

and 

lim  (G(A))‘' =  G.  (28c) 

*  —  r 

To  derive  the  numerically  stable  algorithm  we  modify  equation  (28a)  as 

AT-\n[G{k  -F  l))  =  /ir-ln(G(A)]~ln{exp[/ir-ln(G(A))l}.  (29) 


Defining 


N{k)  ^  /ir-ln(G(A)) 


(30) 


and  substituting  expression  (30)  into  equation  (29)  and  equations  (28)  yields 


A^(A  -F  1 )  =  N(k)  -  ln{exp(/V(A))],  .V(0)  =  A  Tiq 


(31a) 
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G{k  +  \)  =  G{k)txp(N{k)),  G{0)  =  /„  (31b) 

Urn  iV(A:)  =  0,  (31c) 

*  -*  X 

lim  ((/(/:))'' =  exp(/4r)  =  G.  (31d) 

*  ^  X 


Since 


exr'V)^/„  +  3r^[/„  +  U][/„-U]-'~  -,  (32) 

where  A'  is  a  matrix.  The  algorithm  in  equations  (31)  can  be  represented  by  various  forms  as 
follows; 


iV(/t  +  1)  =  tV(A:)-ln[/„  +  /V(>t)],  N{Q)  =  ATq  (33a) 

G(k  +  \)  =  G(k)[I„^N(k)].  G(0)  =  /,  (33b) 

limAf(A)  =  0„  (33c) 

k  —  X; 

lim  (G(k  ))^  =  exp(/t  T)  =  G.  (33d) 

*  -*  ® 

Also, 

N(k  +  \)^  N(k)-H[I„  +  \N{k)][I„-^_N{k)]-'}.  NiO)  =  ATIq  (34a) 

G(fc  +  l)  =  G(A:)[/„  +  iA(A:)][/„-iA(^)]-',  G(0)  = /,  (34b) 

limiV„<)  =  0„  (34c) 

*  -*  X 

lim  (G{k)y  =  exp(/17’)  =  G  (34d) 

*  -•  X 

and 

N{k  +  \)  =  N{k)-\n{[I„  +  \Nik)  +  iiN{kW„-\Nik)  +  iiN{kf]-'},  N(0)  =  ATIq  (35a) 

G(k +  \)  =  G{k)[I,  +  \N{k)  +  i-,N{km,-\N(k)+  ^A(A:)-r',  G(0)  =  /„  (35b) 

\imN(k)  =  0„  (35c) 

k  X. 

lim  {Gik)y  =  exp(^r)  =  G.  (35d) 

k  —  X 


The  matrix  logarithmic  functions  in  equations  (33a),  (34a)  and  (35a)  can  be  computed  using  the 
techniques  shown  in  Section  2.  The  algorithms  in  equations  (33),  (34)  and  (35)  are  numerically 
stable.  The  stability  is  proved  in  Appendix  B. 

Substituting  the  obtained  matrix  G  and  the  available  matrix  A  to  equation  (3c)  gives  the  input 
matrix  H.  Thus,  we  obtain  the  discrete-time  model  in  equation  (3a). 


5  ILLUSTRATIVE  EXAMPLE 

A  practical  example,  which  consists  of  a  stiff  system  matrix,  is  presented  for  finding  the  matrix 
A(G)  from  the  matrix  G(A)  as  follows. 

Consider  the  stiff  continuous-time  system  matrix  of  a  linearized  two  shaft  gas-turbine  model  [19]: 


A  = 


-1.268 

1.002 

0 

0 


-0.04528  1.498  951.5 

-  1.957  8.52  1240 

0  -10  0 
0  0-100 


where  (T (/I)  =  {-1.3417,  -1.8833,  -10,  -100}  and  (t(/1  )r  =  { -0.05367,  -0.07533,  -0.4,  -4} 
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having  a  fixed  sampling  period  T  =  0.04s.  Note  that  \(t{AT)\  >0.5.  The  associated  discrete-time 
system  matrix  with  the  same  fixed  sampling  period  T  =  0.04  s  is 


G  = 


0.9505105993 

0.0375771817 

0 

0 


-0.0016980986 

0.0478132051 

8.967804795 

-0.9246715991 

0.2704783066 

11.7361807569 

0 

0.6703200460 

0 

0 

0 

U.ui83 156389 

where  cr(G)  =  {0.94774510,  0.92743710,  0.67032005,  0.01831564}. 

It  is  desired  to  find  the  matrix  A  from  the  matrix  G  via  the  algorithms  in  equations  (16),  (20), 
(23),  (24)  and  (25).  Since  Re[(7(G)]  >0,  the  square-root  factor  q  can  be  chosen  as  1.  The  infinite 
series  approximation  method  shown  in  equation  (4)  with  y  =31  is  employed  for  finding  the 
approximants  of  the  matrix  exponential  functions  in  equations  (16),  (20),  (23),  (24)  and  (25).  The 
simulation  results  are  shown  in  Fig.  2.  The  algorithms  in  equations  (16)  and  (20)  converge  at  the 
9th  iteration  with  1 1  (^4  (9)  —  )y4  ' '  1 1  =  3.4  x  10'*  and  at  the  5th  iteration  with  1 1  (/I  (5)  —  ^  )/l  " '  1 1  = 

7  X  10"",  respectively,  then  both  algorithms  diverge.  Therefore,  the  algorithms  in  equations  (16) 
and  (20)  are  numerically  unstable.  The  algorithms  in  equations  (23),  (24)  and  (25)  converge  at  the 
9th  iteration  with  \  \{A(9)  —  A)A~'\\  =  4.4  x  10"'®,  at  the  6th  iteration  with  \  \{A(,6)  —  A)A~'\\  = 
2.2  X  10 and  at  the  5th  iteration  with  ||(/1(5)  — /4)/l  "'||  =  3.6  x  10"''',  respectively.  All  three 
algorithms  are  numerically  stable.  Since  the  computational  loads  of  the  proposed  rapidly 
convergent  algorithms  in  equations  (24)  and  (25)  are  greater  than  that  of  the  existing  algorithms 
in  equations  (16),  it  might  be  interesting  to  compare  the  CPU  times  of  the  algorithms  in  equations 
(23),  (24),  (25)  and  (16).  With  the  same  relative  error  level  (1  x  10"^),  the  CPU  times  for  the 
algorithm  in  equations  (23)  with  7  iterations,  the  algorithm  in  equations  (24)  with  5  iterations,  the 
algorithm  in  equations  (25)  with  4  iterations  and  the  algorithm  in  equations  (16)  with  9  iterations 
using  the  VAX  1 1/750  computer  are  0.007,  0.005,  0.28  and  0.14  ms.  respectively.  It  is  observed  that 
the  proposed  algorithms  in  equations  (23)  and  (24)  give  faster  computation  time  than  the  existing 
algorithm  in  equations  (16)  which  involves  two  matrix  inversions  for  computing  the  initial 
condition  and  G"',  while  the  proposed  algorithm  in  equations  (25)  provides  slower  computation 
time  than  the  algorithms  in  equations  (23),  (24)  and  (16)  because  it  involves  four  matrix  inversions 
(i.e.  one  matrix  inversion  for  each  iteration).  From  above  results  we  conclude  that  a  rapidly 
convergent  algorithm  may  not  provide  faster  computation  time  due  to  computational  complexity 
of  the  algorithm.  However,  the  computational  time  due  to  computational  complexity  can  be 
improved  by  using  advanced  computer  technologies  such  as  the  parallel  computational  techniques 
[20],  etc. 


Fi?.  ■  Finding  A  from  G:  O — algorithm  in  equations  (16);  A — algorithm  in  equation  (20); 
□  -ithm  in  equations  (23);  ■ — algorithms  in  equations  (24);  A — algonthms  in  equations  (25) 
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Fig.  3.  Finding  G  from  A: 


■ — algorithm  in  equations  (33);  O — algorithm  in  equations 
□ — algorithm  in  equations  (35). 


(34); 


To  satisfy  the  condition  \a{AT)\  <0.5  in  the  algorithm  in  equation  (16),  we  select  T  =  0.005  s 
and  carry  out  the  simulations  using  algorithms  (16),  (23)  and  (24).  The  algorithms  in  equations 
(16),  (23)  and  (24)  converge  at  the  5th  iteration  with  a  relative  error  2.1  x  10''^  at  the  5th  iteration 
with  a  relative  error  1.3  x  10”  “  and  at  the  3rd  iteration  with  a  relative  error  6.4  x  10“'^, 
respectively.  The  corresponding  computation  times  are  0.14,  0.005  and  0.003  ms,  respectively.  All 
three  algorithms  are  stable;  however,  the  proposed  algorithms  in  equations  (23)  and  (24)  provide 
faster  computation  times  than  the  algorithm  in  equations  (16)  which  involves  two  matrix 
inversions. 

To  find  the  matrix  G  from  the  matrix  A  we  use  the  stable  algoriiluns  in  equations  (33),  (34)  and 
(35)  as  follows. 

Let  In(A')  a:  2R[I,  _  i  ^  where  R^[X-  I„][X  +  /„]' '  and  <?  =  64.  With  the  same 

relative  error  tolerance  ||(G(fc)  —  G)G  "'ll  2: 1  x  10~’,  the  algorithm  in  equations  (33)  converges  at 
k  =4,  while  the  algorithms  in  equations  (34)  and  (35)  converge  at  k  —2.  The  simulation  results 
are  shown  in  Fig.  3.  For  the  purpose  of  comparison,  the  scaling,  squaring  continued-fraction 
method  in  equations  (5),  which  is  a  non-recursive  algorithm,  has  been  employed  to  compute  the 
approximant  G,  ^  with  i  =  5  and  the  same  value  of  g{  =  64).  The  associated  relative  error  is  2  x  10 
The  corresponding  CPU  times  for  the  algorithms  in  equations  (33),  (34)  and  (5)  for  the  same 
q(=64)  are  0.56,  0.42  and  0.1  ms,  respectively,  while  the  corresponding  relative  errors  for  the 
algorithms  are  1  x  10"’,  1  x  10"’  and  2  x  10"’,  respectively  The  non-recursive  algorithm  in 
equations  (5)  provides  faster  computation  time  than  the  proposed  recursive  algorithms  which 
involve  matrix  inversions;  however,  the  proposed  recursive  algorithms  give  smaller  errors  than  the 
non-recursive  algorithm  for  the  same  q.  The  use  of  a  small  q  will  result  in  a  small  round-off  error 
caused  by  computing  the  matrix  (C,)",  where  the  matrix  G,  is  an  approximant  of  the  cxp(AT/q) 
in  equation  (5a)  or  the  matrix  {G^ik))"  where  the  matrix 

G^(k)=  Urn  (G(/t)) 

in  equation  (28c). 


6.  CONCLUSION 

Numerically  stable  and  rapidly  convergent  recursive  algorithms  for  modelling  the  equivalent 
continuous-time  (discrete-time)  model  from  the  available  discrete-time  (continuous-time)  model 
have  been  presented.  The  non-recursive  algorithms  such  as  the  scaling,  squaring  continued-fraction 
method  and  the  square-root,  scaling,  continued-fraction  method  for  the  models  conversion  have 
been  introduced,  and  their  truncation  errors  have  been  analyzed.  Based  on  the  Newton-Raphson 
method,  new  recursive  algorithms  have  been  developed  for  continuous-time  and  discrete-time 
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model  conversions.  It  is  shown  that  the  proposed  algorithms  are  numerically  stable  and  that  the 
convergence  speed  of  the  recursive  algorithm  for  obtaining  the  continuous-time  (discrete-time) 
system  matrix  from  the  discrete-time  (continuous-time)  system  matrix  depends  upon  the  approxi¬ 
mation  of  the  matrix  logarithmic  (exponential)  function,  whereas  the  accuracy  of  the  recursive 
algorithm  depends  on  the  approximation  of  the  matrix  exponential  (logarithmic)  function.  The 
newly  developed  recursive  algorithms  relax  the  constraints  imposed  on  the  existing  model 
conversion  methods  and,  hence,  can  be  applied  to  the  system  which  may  contain  both  large  and 
small  eigenvalues  with  a  fixed  sampling  period.  A  practical  example  has  been  presented  to 
demonstrate  the  effectiveness  of  the  proposed  procedures.  Thus,  utilizing  the  proposed  model 
conversion  techniques,  we  can  analyze  and  design  a  sampled-data  system  via  the  well-developed 
continuous-time  or  discrete-time  techniques  [2,21], 

It  should  be  noted  that  the  proposed  algorithms  do  converge  rapidly;  however,  their  com¬ 
putational  complexity  is  greater  than  the  existing  model  conversion  algorithms.  As  a  result,  the 
proposed  algorithms  may  not  necessarily  give  the  fast  computational  time.  The  computational  time 
often  depends  upon  the  convergence  speed  and  the  computational  complexity  of  the  algorithm.  The 
computational  time  due  to  the  computational  complexity  can  be  improved  by  using  advanced 
computer  technologies  such  as  the  parallel  computational  techniques  etc. 
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APPENDIX  A 

Proof  of  the  Numerical  Stability  of  the  Algorithm  in  Equation.^'  (21) 
Consider  the  algorithm  in  equations  (23)  with  q  =  1  and  M{k)  exp(- A(k)T)G,  i.e, 

M{k  +  l)  =  M(k)np[l,-M(k)].  M(0)  =  G 


(A.  la) 
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Aik  +  l)  =  Aik)-j.il,-Mik)].  /t(0)  =  0. 

(A.  lb) 

lim  Mik)  =  /, 

(A.lc) 

k  —  ■Xi 

lim  Aik)  =  A. 

(A.ld) 

it  is  desired  to  show  that  the  iteration  in  equations  (A.I)  is  numerically  stable  in  the  sense  that  a  small  perturbation 
introduced  at  the  >tth  iteration  due  to  round-off  errors  cannot  lead  to  unbounded  perturbations  on  succeeding  iterates. 
Let  the  perturbed  models  be  M(k)  and  ^[k)  and  the  associated  round-off  errors  be  E(k)  and  F(k),  respectively.  Hence 


M(k)  =  M(k)  +  E{k)  (A.2a) 

A{k)T  =  A(k)T  +  F(k).  (A.2b) 

Our  aim  is  to  analyze  how  the  error  matrices  E{k)  and  F{k)  propagate  at  the  (k  -t-  1)  stage.  To  simplify  the  analysis  we 
assume  that  no  round-off  errors  occur  when  we  compute  M{k  +  1)  and  -l-  1)  in  the  following: 

W(L I)  =  A?(*)exp[/, -/!?(*)]  (A. 3a) 

Mk  +  i)T  =  A(k)T M(k)].  (A.3b) 


Substituting  equations  (A. 2a)  into  (A. 3a),  expanding  the  matrix  exponential  function  in  equation  (A. 3a)  and  omitting  the 
high  order  trivial  terms  of  E[k)  we  have 

£(/t  -f  1 )  =  A?(/t  I )  -  A/(L  -(-  I )  =  [M(k)  +  E(k )]{/.-(-  [/,  -M(k)-  E{k )] 

-(-  I  [(/,  -  M(k)f  -  (/,  -  M(k))E(k)  -  E{k)(I,  -  A/(A:))1  -I-  •  •  ■  i  -  /V/(*)  exp[/,  -  M{k)] 

=  M(k)[E  +  (/.  -  M(k))  +  'i(I.  -  M{k)f  +  ■  ■  ■] 

-t-  E(k)[l.  +  (/,  -  m ))  + 1  (/, -  M(k ))'  4-  •  •  •  ] 

-  M{k)[E(k)  +  i  [(/,  -  M{k  ))E(k )  +  E(k){I,  -  M(,k))]  +  ■  ■  ■]  -  M{k)  exp(/,  -  M(k )[ 


=  £(*)exp{/,  -  M(k)]  -  M(k){E(k)  +  ^[{1,  -  M(k))E{k)  +  E(k){I,  -M(k))]  +  ...}  (A.4a) 

Also,  from  equations  (A. lb),  (A.2b)  and  (A.3b)  we  obtain 

F(k  +  1)  =  A(k  +  1)7  -  Aik  H- 1)7  =  Aik)T  -l,+  !^(k)  -  A{k)T  +  1,~  Mik)  =  f(lfc)  +  £(L).  (A.4b) 

When  k  -*  <X),  exp(  — ,4(A:)7) -•  Cf"',  hence  M(k)  =  exp(  —  A(k)T)G  As  a  result,  exp[/,  -  A/(/t)) Thus,  equa¬ 
tions  (A.4)  becomes 

£(* -t-l)  =  0,  (A.5a) 

Fik  +  i)  =  F(k)  +  Eik).  (A.5b) 


The  discrete-time  state  equation  in  equation  (A.Sb)  with  an  identity  system  matrix  is  stable.  If  we  make  a  further 
assumption  that  no  round-off  errors  are  introduced  at  the  ik  +  2)  stage  of  the  iteration,  then  equation  (A.Sb)  becomes 
F(k  4-  2)  =  Fik  4-  1)  4-  £(*:  4- 1)  =  Fik  +  1).  This  suggests  that  the  perturbation  introduced  at  the  Ath  stage  has  only  a 
bounded  effect  on  succeeding  iterates.  Therefore,  the  algorithm  in  equations  (A.  1)  is  numerically  stable.  In  a  similar  manner, 
we  can  show  that  the  algorithms  in  equations  (24),  (25)  and  (26)  with  q  >  I  are  numerically  stable.  Note  that 
Mik)  =  exp(  — /4(lc)7TC  i;  G*'G  =  even  if  the  given  matrix  G  is  a  stiff  matrix. 
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Proof  of  the  Numerical  Stability  of  the  Algorithm  in  Equations  (31) 
Consider  the  algorithm  in  equations  (33)  with  q  =  I  and  Nik)  ^  AT  —  ln(G(it)),  i.e. 


Nik  +  l)  =  Nik)-ln[f  +  Nik)].  NiO)  =  AT  (B.la) 

Gik  4-  1)  =  Gik)[f  4-  Nik)].  G(0)  =  /,  (B.lb) 

lim/V(lc)  =  0,  (B.lc) 

lim  G(<: )  =  exp(/<  7)  =  G.  (B.ld) 

*  -•  dC 


It  is  desired  to  show  that  the  iteration  in  equations  (B.l)  is  numerically  stable  in  the  sense  that  a  small  perturbation  at 
the  Ath  iterate  due  to  round-off  errors  cannot  lead  to  unbounded  perturbations  on  succeeding  iterates. 

Let  the  perturbed  models  be  ^(k)  and  Gik)  and  the  associated  round-off  errors  be  Eik)  and  Fik),  respectively.  Hence 

i}(k)  =  Nik)  +  Eik)  (B.2a) 

G(k)  =  Gik)+  Fik).  (B.2b) 


Our  purpose  is  to  analyze  how  the  error  matrices  E(k)  and  F(k)  propagate  at  the  (*  +  l)th  stage.  To  simplify  the  analysis 
we  assume  that  no  new  round-off  errors  occur  when  we  compute  Nik  +  1)  and  Gik  +  1)  in  the  following: 

!iikA-\)  =  f}ik)-\n[[,  +  f}ik]]  (B.3a) 

Gik  +  \)  =  G(k)[I,+  f}ik)].  (B.3b) 

Substituting  equations  (B.2)  into  (B.3),  expanding  the  matrix  logarithmic  function  in  equation  (B.3a)  and  omitting  the  high 
order  trivial  terms  of  Eik)  yield 

Eik +  l)  =  E(k)- {Eik) -'i  [Nik )£(* )  4-  Eik )Nik )]  4- 1  [Nik  fEik )  4-  Nik  )Eik )iV(*: )  4-  £(*  )Nik )2) - }  (B.4a) 
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and 

F{k  +  \)  =  F(k)  +  G(k)E(k)  +  F(k)N(k).  (B,4b) 

When  k  -»  oo,  G(k)  exp{AT),  hence  N{k)  =  AT  —  ln(G(/c))  =  AT  -  ln(exp(/< D)  -•  0,.  Thus,  equations  (8.4)  becomes 

£(*  +  l)  =  0,  (B.Sa) 

£(*  + !)  =  /■(*)  +  <?(*  )£•(*).  (B.5b) 

The  discrete-time  state  equation  in  (B.5b)  with  an  identity  system  matrix  is  stable,  ir  we  make  a  further  assumption  that 
no  new  round-off  errors  are  introduced  at  the  (k  -i-  2)  stage  of  the  iteration,  then  equation  (B.5b)  becomes 
F(k  -I-  2)  =  F(k  -b  1)  -(-  G(^  +  -b  I)  =  F(k  +  1).  This  suggests  that  the  perturbations  introduced  at  the  k  th  stage  have 
only  a  bounded  effect  on  succeeding  iterates.  Thus,  the  algorithm  in  equations  (B.l)  is  numerically  stable.  In  a  similar 
manner,  we  can  prove  that  the  algorithms  in  equations  (32),  (33)  anJ  t34)  with  q  >  1  are  numerically  stable.  Since 
N(k)  =  AT  —  ln(G{A:))  s:  AT  —  AT  =  0„  the  given  matrix  A  can  be  a  stiff  matrix. 


